Unjamming and emergent nonreciprocity in active ploughing through a compressible viscoelastic fluid

A dilute suspension of active Brownian particles in a dense compressible viscoelastic fluid, forms a natural setting to study the emergence of nonreciprocity during a dynamical phase transition. At these densities, the transport of active particles is strongly influenced by the passive medium and shows a dynamical jamming transition as a function of activity and medium density. In the process, the compressible medium is actively churned up – for low activity, the active particle gets self-trapped in a cavity of its own making, while for large activity, the active particle ploughs through the medium, either accompanied by a moving anisotropic wake, or leaving a porous trail. A hydrodynamic approach makes it evident that the active particle generates a long-range density wake which breaks fore-aft symmetry, consistent with the simulations. Accounting for the back-reaction of the compressible medium leads to (i) dynamical jamming of the active particle, and (ii) a dynamical non-reciprocal attraction between two active particles moving along the same direction, with the trailing particle catching up with the leading one in finite time. We emphasize that these nonreciprocal effects appear only when the active particles are moving and so manifest in the vicinity of the jamming-unjamming transition.


Supplementary Note 1. AGENT-BASED SIMULATIONS
To ensure that the assembly of soft particles of area fraction φ remains in a positionally disordered state, we work with a modified 2D Kob-Andersen 65 : 35 (A:B) binary mixture [1,2], a prototype glass-forming liquid. Particles i and j, interact via a 2-body soft repulsive potential, where r ij is the distance between the particles i and j. Our choice of parameter values, v 0 = −112 ij , v 2 = 192 ij , v 4 = −84 ij and the cut-off r c,ij = σ ij , ensures that both the potential and the force is smooth at the cut-off distance and the interaction is purely repulsive. We fix the energy and length scales to be in the units of AA and σ AA , respectively by setting AA = 1.0, σ AA = 1.0. The other parameters are chosen following the Kob-Andersen model [1] i.e. BB = 0.50, σ BB = 0.88, AB = 1.50, σ AB = 0.80.
Of these a small fraction of particles φ a is made active -their dynamics is described by active Brownian particles (ABP) [3][4][5][6][7][8] immersed in a background of passive particles. All particles are subject to a thermal noise ϑ of zero mean and variance equal to 2γT (setting k B = 1), obeying FDT. The subset i ∈ A of ABPs are subject to additional active stochastic forces f i = f n i ≡ f (cos θ i , sin θ i ). The orientation of the propulsion force θ i undergoes rotational diffusion, described by an athermal noise ξ i , with zero mean and correlation ξ i (t)ξ j (t ) = 2τ −1 δ ij δ(t − t ). Its effect on the particle-dynamics appears as an exponentially correlated vectorial noise with correlation time τ , which being unrelated to the drag γ, is not constrained by fluctuation-dissipation relation.
The full dynamics is described by the Langevin equation, where 1 (i∈A) is the indicator function which ensures that the active forces are only imposed on particles i belonging to the active set A.
To proceed with the numerical simulation, we take the overdamped limit (drop inertia) and convert the resulting dynamical equation to non-dimensionless form by scaling, and setting σ AA = AA = 1, to arrive at, Similarly, the orientational dynamics (which is purely diffusive), can be written as,θ We perform Brownian dynamics (BD) simulations at fixed particle-number, volume of the system and temperature of the heat-bath (NVT) in 2-dimensions using a square box of reduced length L 0 = 45 and 90, with periodic boundary conditions (PBC). We use predictor-corrector algorithm (Euler-Trapezoidal method) and forward Euler method to update the positions and orientations, respectively [9,10]. For all the simulations related to Figs. 1, 2 (main text) we keep the area fraction of active particles fixed at φ a = 0.017 (dilute limit), and vary the number of passive particles constituting the medium to control the overall density or area fraction φ. The simulations to generate data for Figs. 3, 4, 5, 6 (main text) are performed with one active particle in the passive medium. For different values of the area fraction φ, we first equilibrate the passive medium at the reduced temperature T = 0.5. We then monitor the approach of the passive medium towards a glass transition, by computing the relaxation of its density fluctuations, measured by the two point overlap correlation function ( Fig. Supplementary Figure 1A); the transport of tagged particles, measured by the mean squared displacement (Fig. Supplementary Figure 1C) and its rheological properties, as we systematically vary φ. Only in this section, τ indicates the time elapsed starting from the time-origin, t, and is not to be confused with the τ used in all the other sections of this Supplementary Information as well as in the main text, representing the persistence time of the active particles.
Overlap function, Q(τ ), a measure of the density relaxation, is defined as with ω(x) = 1 for x < 0.3 and zero otherwise. The ... sign indicates averaging over different time-origins (t), and different realisations of the system. The α-relaxation time, τ α is calculated by setting Then τ α , for different area fractions φ, is fit to the Vogel-Fulcher-Tammann (VFT) form ( Fig. Supplementary Figure 1B), to obtain φ VFT , the VFT glass transition area fraction. Mean squared displacement (MSD), is defined as, The . . . denotes averaging over the time-origins and multiple statistically independent realizations. The long time diffusion coefficient (D ∞ ) is computed from the long time limit of the MSD, Passive microrheology, is a measure of the viscoelastic(shear) properties of the medium, and can be obtained from the frequency transformed MSD, ∆r 2 (ω) [11], where, Here, 'a' is the radius of the particle and Γ(x) is the Γ-function. We choose T = 0.5 for all the simulations related to Fig. 1 as (1) at this temperature cage-hopping dynamics is clearly visible for higher area fractions (i.e., φ = 0.7−0.8) of the medium, and (2) for moderate area fractions (i.e., φ = 0.4 − 0.6) the temperature is not too close to the glass transition temperature where the dynamics becomes extremely slow and it becomes very difficult to acquire enough data within a reasonable run-time. Both the equilibration and data-generation runs are for a time ≥ 100 τ α .

A. Creation of a density map from agent-based simulations
We use the instantaneous configuration (i.e., position-coordinates of all the passive particles) to generate a coarsegrained density field that has been used in Figs. 3-5 in the main text. For this, we divide the entire simulation box of size L 0 into a square grid of grid length L g , and number of cells N g = ( L0 Lg ) 2 . For a typical cell with indices (k, l), we compute the density field using the following formulae where d = |x i − x k,l | is the distance of the i-th particle from the centre of the grid (x k,l ) and d 0 is the coarse-graining length scale, and N is the normalisation constant. We have used different values for L g and d 0 , ensuring that d 0 > L g and that both are ∼ O(1). For Fig. 3, 4 since we compare the density fields between different area fractions, we choose which ensures that (S14) gives us the normalised probability density of the passive medium. For Fig. 5 and Fig. 6 (main text), we choose which ensures that (S14) gives us normalised number density.

S6
B. Identification of under-dense regions and analysis of their geometry We use the density profile that we discussed in Sect. Supplementary Note 3 A and identify the low density regions by employing an upper cut-off (∼ 4 × 10 −4 ) on the grid defined above. The results described in Fig. 4A in main text are robust to reasonable variation of this threshold. Any grid point (cell) that has density less than or equal to the aforementioned threshold is considered as low density region. This low density regions are represented by a set of points defined on the grid z i , where i ∈ [1, N u ] where N u is the total number of under-dense grid points. We then compute moment of gyration tensor,

Supplementary Note 4. LINEARISED HYDRODYNAMIC THEORY
We first write the full dynamics of the medium and the active particles, and their interplay. The medium is described by a local density field ρ(r, t) and a velocity field v(r, t), and obeys, the continuity equation for ρ ∂ t ρ + ∇ · (ρv) = 0 (S17) and force-balance equation for v, In (S18) above, B > 0, since the pressure due to the medium opposes the local forcing from the active particles. The dynamics of the active particles in the dilute limit is given by where γ is the drag on the active particle and C is the extent of back reaction of medium on the translational dynamics of the active particle. n i ≡ (cos θ i , sin θ i ), and the athermal orientational noise ξ i has zero mean and is delta-correlated, In general, the frictions Γ, γ, the active force magnitude f and coefficients B and C could depend on ρ.
A. Single active particle moving through the compressible medium in d = 2 We will now restrict ourselves to a single active particle moving through the compressible medium, where R(t) is the coordinate of the active particle. As the active particle moves through the medium, it creates density inhomogeneities, which relax over time; we will be interested in a linear analysis about the uniform density of the ploughed medium. In this limit, we take Γ and B to be independent of ρ, and ignore the back-reaction term C. Plugging (S21) into the continuity equation (S17), we get where we have absorbed Γ and redefined B , f as B/Γ and f /γ respectively. To simplify analysis, we will assume that the persistence time of the active particle is large compared to the timescale of observation, i.e., n(t) is independent of time. After initial transients, the density of the medium can then be written in terms of the collective coordinate, Transforming coordinates to the co-moving frame of the active particle (without loss of generality, we take n to point along thex-direction), and so, where v 0 ≡ f /γ. The equation for the transformed density ρ(u, w) in the co-moving frame now reads, where ∇ ≡ (∂ u , ∂ w ).
We now express the local density as a small deviation (positive or negative) from the uniform ρ 0 , and rewrite the above equation to linear order in this small deviation, i.e., we set ρ → ρ 0 + ρ in (S25), and linearise to get, This linear equation can be solved by fourier transforming (u, w) → (q u , q w ), and then evaluating the inverse fourier transform using contour integration. To proceed, we first note that, which leads to, where q 2 = q 2 u + q 2 w , from which we get,

S8
We now obtain ρ(u, w), by inverse fourier transforming ρ(q u , q w ), where I is an integral over q u , keeping q w fixed, We evaluate I in (S32) using contour integration. For this, we determine the simple poles q ± of the integrand by factorization, where the length-scale ξ is given by, when we restore the original parameter definitions. Note that the sgn Im(q + ) = +1 and sgn Im(q − ) = −1. Thus we close the contour of integration in the upper half plane around q + when u > 0, and in the lower half plane around q − when u < 0. FIG. Supplementary Figure 3. Choice of contours for A u > 0 and B u < 0.

B A
Using, q + − q − = 2i ξ −2 + q 2 w , we can evaluate the residues, and The integral I (S32) can now be read out. For u > 0, For u < 0, we can write the above expressions in compact notation where u ± denotes positive and negative u, respectively. Note that the function in (S41) is even in q w . We now evaluate the integral in (S31), cos (wq w ) dq w = I 2 (u, w) for u, w > 0 where K 0 and K 1 are the modified Bessel functions of the second kind.
Therefore, for the whole range of u, where, The fore-aft asymmetry of the density profile is most apparent when we set w = 0 + . In this case, for u > 0, Now using the asymptotic expansions [13] and for large z, we see immediately that the excess density profile in the moving frame is fore-aft asymmetric. While the profile decays exponentially following a pile up in the region u > 0, there is a slower power-law decay in the region u < 0, going as |u| −3/2 .
In a similar way, we obtain the density profile in the direction transverse to the direction of motion, B. Two active particles moving through the compressible medium in d = 2 We will now be interested in how two active particles at R 1 , R 2 , moving through the compressible medium, affect each others dynamics. Let us again, take the limit of large persistence time -in this limit the orientations n 1 and n 2 are time independent. To O(ρ), these particles, individually, leave a time dependent asymmetric wake described by ρ(u, w, t) in (S46) To compute the back reaction of this wake on the movement of the active particles, we need to include terms to O(ρ) in (S20). Thus, These equations may be cast in terms of the relative coordinate R rel = R 1 − R 2 and the centre of mass R cm = (R 1 + R 2 )/2. We look at the following scattering geometries and initial conditions -1. Active particles 1 and 2 oriented along the fixedx direction (n 1 = n 2 =x), with initial positions along the x-axis (particle 1 ahead of particle 2) Note that the two active particles are identical, and so their self contributions to this order are the same. Therefore, in terms of the relative coordinate X rel = X 1 − X 2 > 0 and the centre of mass X cm = (X 1 + X 2 )/2, (S53), (S54) can be recast as, and We focus on the relative coordinate X rel . The effect of particle 2 on 1 can be obtained from (S48) The effect of particle 1 on 2 is obtained from (S49), where A = f (ρ0+ρ(0,0)) Since the two active particles are identical, ρ 1 = ρ 2 = ρ, and sȯ This dynamical equation has a stable fixed point at X * rel ≈ 1.5064ξ (indicated in the inset of Fig. 7C). This means that, according to this linearised 2-particle theory, there will be a bound 2-particle state, i.e., eventually both particles will move forward keeping the same relative distance.
In the limit z → ∞, the equation for the dynamics of X rel becomeṡ Cρ 0 A γξ Similarly, for the motion of the centre of mass, we have, Cρ 0 A γξ (S60) together with (S59) implies that if the two particles start from an initial separation X rel ξ, the separation decreases with time with an initially increasing speed of approach, which, over time, decreases and the particles, eventually, move with the same speed, being bound to each other at a finite separation, X * rel . This is a signature of nonreciprocal sensing.
We may extract the scaling behaviour from the dominant term in (S60). When X rel ξ, one may obtain, as stated in the main text.
2. Active particles 1 and 2 oriented along the fixedx direction, with initial positions along the y-axis (particle 1 to the left of particle 2) and The symmetry of the problem ensures that the last term in (S64) will be zero. To evaluate Y rel , we compute and for Y rel ξ → 0, Unlike in Case (1), particles far separated along y do not feel a dynamical attractive force. However, if the y-distance between particles is of order ξ, then they will feel an effective attraction, which will make their trajectories converge towards each other. At late times, when Y rel ξ, one may obtain the scaling behaviour The linearised theory shows that the excess density in d = 1 is very different from higher dimensions, because of the absence of a back flow in the medium. We will verify this against exact numerical solutions of the nonlinear equations in d = 1.
We first recall the linearised equations written in the comoving coordinate u of the motile particle, As before, the solution is obtained by Fourier transformation, (S71) S14 which upon inverse fourier transforming gives, The integrand has a simple pole at where the length scale ξ = Bρ 2 0 /(v 0 Γ). Note that the sgn Im(q + ) = +1. For u > 0, the integrand will tend to zero only when Im(q u ) > 0. For u < 0, the integrand will tend to zero only when Im(q u ) < 0. Thus, the contour of the integral should traverse the upper half plane around q + when u > 0, and in the lower half plane when u < 0.
This immediately shows that the excess density Unlike the 2d density profile, there is no trailing wake in d = 1. We numerically solve the 1-dimensional non-linear equations corresponding to (S17)-(S20), with no noise. We have used finite volume discretization with an exponential scheme [14] for the convective flux as implemented in NIST-FiPy [15]. For the stability of the numerical scheme, we consider a diffusion term in the density field evolution D∂ 2 x ρ with a very small diffusion constant (D ∼ 10 −3 ). According to the exponential scheme, the convective problem is discretized and solved in each grid using an approximate form of the convective transport with constant coefficients ∂ x (vρ − D∂ x ρ) = ∂ x J x = 0. Thus, the local (over a grid) analytical solution for the density becomes exponential. This analytical expression is then used to estimate the flux. The flux at the interface of the i th and (i − 1) th grid in this scheme is given by where P e = v∆x D is the Péclet number and ∆x is the grid spacing. This scheme guarantees positive solutions and has low diffusive error as the flux is formulated using the exact solution. A first order temporal discretization was used in combination with a sweep between each iteration using Newton's method. Numerical solution of the nonlinear equation shows perfect agreement with the linear theory for d = 1 and that the density profile has the exponential decay in the front without any trailing wake, consistent with the linear theory ( Fig. Supplementary Figure 5A). We also capture the dynamical transition (Fig. Supplementary Figure 5B) of the active particle with increasing propulsion force f .

Supplementary Note 7. FINITE-SIZE EFFECTS
We observed a strong finite size effect while analyzing the simulation data. This is because there is an intrinsic length scale set by ξ whose effect on the asymptotic decay can be corrupted due to the effects of periodic boundaries if the system size is not sufficiently large. This implies that one should be careful while comparing the theoretical S15 predictions of the linearised hydrodynamic theory with simulations. For instance, at a system size of L 0 = 45σ AA with N = 1260 particles (shown in red in Fig. Supplementary Figure 6), the decay of the density wake shows a (u/ξ) −1/2 scaling (Fig. Supplementary Figure 6B). It is only at a sufficiently large system size L 0 = 90σ AA with N = 5040 particles (shown in blue in Fig. Supplementary Figure 6), that we obtain a (u/ξ) −3/2 scaling (Fig. Supplementary  Figure 6C), that is consistent with the linearised hydrodynamic theory.